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Abstract 

We consider the linear regression model with a spatially corre- 
lated error term on a lattice process. When we estimate coefficients 
in the linear regression model, the generalized least squares estimator 
(GLSE) is used if the covariance structures are known. However, the 
GLSE for large spatial data sets is impractically time-consuming be- 
cause it includes the inversion of the covariance matrix of error terms 
in different spatial points that is the size of the number of observations. 
To reduce the computational complexity, we propose the pseudo best 
estimator (PBE) using spatial covariance structures approximated by 
separable covariance functions. We derive the asymptotic covariance 
matrix of the PBE and compare it with those of the least squares 
estimator (LSE) and the GLSE through some simulations. They also 
imply that the effect of the misspecification of the covariance matrix 
for the GLSE is examined. Monte Carlo simulations demonstrate the 
improvement of the LSE, which does not contain the information of the 
spatial covariance structure, by the PBE using separable covariance 
functions even if the true process has an isotropic Matern covariance 
function. Additionally our proposed PBE is computationally efficient 
relative to the GLSE for large spatial data sets. 
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1 Introduction 



Recently various statistical methods for spatial data has been investigated. 
Among them, a linear regression model with a spatially correlated error term 
has been playing an important role in a wide variety of scientific fields such 
as geostatistics, econometrics and forestry. 

To estimate coefficients in the linear regression model, we often use the 
generalized least squares estimator (GLSE). However, the GLSE for large 
spatial data sets is impractically time-consuming because it includes the in- 
version of the covariance matrix of error terms in different spatial points n 
which requires 0(n 3 ) operations. For example, Isaaks and Srivastava (1989) 
and Huang et al. (2010) introduces the Walker Lake data set and the soil 
moisture index derived from other GIS layers which are observed on a spatial 
lattice of 260 x 300 and 100 x 100 regular grid points respectively. There is 
a rich literature on the analysis of large spatial data sets (Furrer et al. 2006; 
Banerjee et al. 2008; Cressie and Johannesson 2008; Lindgren et al. 2011). 

To deal with this problem, we propose the pseudo best estimator (PBE) 
using spatial covariance structures approximated by the covariance function 
which is expressed by the product of the covariance functions of the causal 
autoregressive process. The covariance matrix of the PBE has the separable 
structure and it facilitates computational procedures for large spatial data 
sets by the property of the Kronecker product. Genton (2007) applied the 
separable approximation to the prediction of large space-time data sets by 
solving the nearest Kronecker product for a space-time covariance matrix 
problem and obtained the good predictive performance. 

In this paper, we derive the asymptotic covariance matrix of the PBE for a 
lattice process and compare it with those of the least squares estimator (LSE) 
and the GLSE in Yajima and Matsuda (2008) through some simulations. The 
LSE is considerable efficient in terms of the fast calculation of the estimator 
for large spatial data sets and Yajima and Matsuda (2008) obtained the 
necessary and sufficient conditions where the asymptotic covariance matrix 
of the LSE is identical with that of the GLSE. However, our simulations 
illustrate that the PBE improves the accuracy of the LSE substantially when 
these conditions are not satisfied and shows a good performance as well as 
the LSE even if these conditions hold. In simulations, the difference in the 
mean squared error between the PBE approximated by separable covariance 
functions and the GLSE is small even in the case of the true process with 
the isotropic Matern covariance function. Additionally, our proposed PBE is 
computationally efficient relative to the GLSE. 

Our work can be regarded as an extension of Amemiya (1973) and Engle 
(1974) which considered the asymptotic properties of the GLSE when the 
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covariance structure of the true process was incorrectly identified in time 
series literature. The GLSE is infeasible unlike the PBE because the true 
covariance matrix is unknown in practice. Our paper also examines the effect 
of the misspecification of the covariance function for the GLSE in spatial 
statistics. Koreisha and Fang (2001) investigated the finite accuracy of the 
GLSE and the PBE in time series case. 

The subsequent sections are organized as follows. We introduce a linear 
regression model with a spatially correlated error term and propose the PBE 
by separable covariance functions in Section 2. In Section 3, we derive the 
asymptotic covariance matrix of the PBE and introduce the necessary and 
sufficient conditions where the asymptotic covariance matrix of the LSE is 
identical with that of the GLSE. In Section 4, computer experiments are 
conducted to see the finite sample performance of our asymptotic result and 
compare the asymptotic variance of the PBE with that of the LSE. A con- 
clusion and future studies are mentioned in Section 5. Technical proofs of 
the lemmas and the theorem are given in Appendices A and B. 

2 Linear regression model and some estima- 
tors 

From now on for simplicity we consider the sampling region of the square on 
the plane. Define Z = {0, ±1, ±2, . . .}. Let Z 2 be the integer lattice points 
in the 2-dimensional Euclidean space. For t = (ti,^) (£ Z 2 ), consider the 
regression model of the form 

Vt = x' t p + e t , 

where {yt} is an observed sequence, X t = (xt,i, ■ ■ ■ , xt,p)' is a p-vector of 
nonstochastic regressors, (3 — (Pi, ... , j3 p ) is a vector of unknown regression 
coefficients and the prime denotes the transposition. Hereafter it is assumed 
that (y t , X t ) is observed on the sampling domain P N = {t = (t 1: t 2 )' £ Z 2 |l < 
ti < N, 1 < t 2 < N}. The error terms {et} follow a stationary random field 
with mean and the spectral density function /(A), A = (Ai, A2) G [— 7r, ti] 2 . 
Then the covariance function %(h) of {e t } is given by 

7 e (fc) = E[e t e t+h \ = / exp(iti \)f(\)d\, 
Jii 2 

where II = (— ir, ti], h = (hi, h 2 )' (E Z 2 ) and h X = hiXi + /i 2 A 2 . 

When we estimate coefficients in the linear regression model, the GLSE 
is given by 
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where X = (X"(i,i), ■ • • , ^(i,n), -^(2,1), ■ ■ ■ , X( 2 ,n), ■ ■ ■ , X(n,i), ■ ■ ■ , ^(n,n)) , 
V = (2/(1,1) , ■ • • , V(i,N), 2/(2,1), • • • , V(2,N), V(N,i), U(n,n)) , E = £[ee'] and 

e = ( e (l,l)) • • • 5 e (l,N), e (2,l), • • • , e (2,AT), • • • , C(JV,1), • • • , e (N,N)) • 

It is known that the GLSE is the best linear unbiased estimator (BLUE). 
However, if the true covariance structure is unknown, the GLSE is infea- 
sible. Additionally, the operation count for computing E -1 (N 2 x N 2 ) of 

Pglse (= {X'H^Xy 1 X'T,- 1 y^j is of order N 6 . Hence as the sample size 
is larger, the computation becomes a more formidable one in practice. For 
example, the Walker Lake data set in Isaaks and Srivastava (1989) consists 
of two variables measured at 78000 points on a spatial lattice of 260 x 300 
regular grid points. In Huang et al. (2010), there is the soil moisture index 
derived from other GIS layers on 100 x 100 points. 

To reduce the computational burden, we consider the approximation of 
the true covariance function by the product of the covariance functions of the 
causal autoregressive process of order P (AR(P)) in time series literature. 
Then we obtain the following estimator 

Ppbe= (x't^Xy 1 X't- l y, 

where S = Ei S 2 and Ej is the covariance matrix of causal AR(Pj) (i = 
1,2). This kind of estimator, in which the true covariance matrix is replaced 
with the incorrect one, is called the pseudo best estimator. Each element 
of E is denoted by the separable covariance function 7(^1,^2) = 7i(^i) x 
72(^2) and '-fi(hi) represents the autocovariance function of AR(Pj) (i = 1,2). 
Note that hi and h 2 of 7(/ii, h 2 ) correspond to those of %(hi, h 2 ). From the 

property of the Kronecker product, E _1 = ^E x E 2 j = Sf 1 E 2 1 and 

we can obtain the exact form of the inverse matrix of the covariance one Ej 
(i = 1,2) of the autoregressive process (e.g., Anderson 1974; page 576), so 
that it is much faster to calculate ftp be by the separable approximation of 
the true covariance function. 

The LSE is also the alternative candidate for the GLSE because it does 
not include the inversion of the covariance matrix and is defined by 

P LSE = (X'xy 1 X'y. 

We will compare the accuracy of the GLSE, the PBE and the LSE by eval- 
uating the asymptotic covariance matrix. 
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3 The asymptotic properties of Pglsei Plse 

and $ PBE 

In this section, the asymptotic covariance matrices of Pglse, Plse and (3pbe 
are derived. First we introduce some assumptions. 
Define 

N-hi N-h 2 

t 1= i t 2 =i 

N N 

— X {ti+hi,t2+h 2 ),i x (ti,t 2 ),j, hi,h 2 <0, 

tl=l—hl *2 = 1 — /l2 

AT N-h 2 

= Y ^2 X (ii+hi,t2+h 2 ),i x (t 1 ,t 2 ),j, h 1 <0,h 2 >0, 

t-L=l-hi *2 = 1 

N-hi N 

= Y ^2 x (ti+hi,t2+h 2 ),i x (ti,t2),j, hi>0, h 2 <0. 

tl=l t 2 =l-h 2 

(a) 0^^(0,0) ->■ oo as N ^ oo,i = 

(b) lim JV _ oo af + ' ,1 ' JV+ ' ,2) (0, 0)/af ,JV) (0,0) = 1 for every i and h u h 2 , i = 
1, . . . ,p and /ii, h 2 G Z. 

(c) The limit of 

J N > N )(h h \ 



{a<r>(0,0)a<f>,0)}^ 



llj 

-if/ i i i i i if, 



as iV — > oo exists for every i,j and /ii, /12, (i, j — 1, • • • ,p and /ii, /t,2 £ Z). 
Let 

Pij(h 1 ,h 2 ) = lim i\f' N \h u h 2 ) 

and ^2) be the pxp matrix with («, j)th element Pij(hi, h 2 ). 

(d) i?(0, 0) is nonsingular. 

(e) {e t } is a unilateral MA process, 

00 00 

e {tite) = ^2 Ql,mV(ti-l,t 2 -m), 

1=0 m=0 

where {0/, m } satisfies J2i, m \^i,m\ 2 < 00 with o ,o = 1 and {^( tl ,t 2 )} is a white 
noise with Var(r] itlM) ) = a*. 
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Moreover if we define 9(z 1 ,z 2 ) = J2Zo E™=o ®i,mA z 2 and 0(^1, 2 2 ) = 
0(z 1 ,z 2 )- 1 = EwEm=o^,™^2 m . then Wi,™} satisfies J2 l>m \(f)i, m \ < 00 
with 0o,o = 1- 

(f) /(A) is a positive continuous function in [— tt,ti] 2 . 

(g) N m&xi< tut2 < N (x( tut2)i i) 2 /alf'^ (0,0) ->■ as N ->■ 00 (i = 1, . . . ,p). 

(h) /(Ai, A 2 ) = /(-Ai, A 2 ) = /(Ai, -A 2 ) = /(-Ai, -A 2 ). 

We make brief comments on the assumptions. We can view (a)-(d) as a 
two-dimensional version of so-called Grenander's conditions on X t . Under 
(a)-(d), there exists a Hermitian matrix function M(Ai,A 2 ) with positive 
semidefinite increments such that 

R(h l ,h 2 ) = / exp^OxAi +/i 2 A 2 ))dM(Ai,A 2 ) 
(see Cohen and Francos 2002). Put 

(N N \ / N N \ 

E E ^, t2) , k ^ XM) EE *( tl , t2 )^ (tiAi+feA2) 
t!=it 2 =i J \ti=i< 2 =i y 

/ ((2 I )=(4r , (0,0)a«r ) (0,0)) 1/2 ) 

and define 

/Ai ;>A2 
/ m < j^' N \u 1 ,u 2 )doj 1 duj2. 
-TT J —TT 

Let M k i(\i,X 2 ) and M^' N \\i, A 2 ) be p x p matrices with (/c,Z)th element 
of M(Ai, A 2 ) and M^ N ' N ^(X 1: A 2 ) respectively for /c, / = 1, . . . ,p. If we regard 
M < ^ N ' N \X 1 , A 2 ) and M(Ai,A 2 ) as matrix measures in IT 2 , R^ N \h u h 2 ) = 
[lki' N \hi,h 2 )\ and R(hi,h 2 ) are their characteristic functions respectively. 
Then (c) implies 

M^(A 1) A 2 )4M(A 1) A 2 ) 

as N — > 00, where means M^ N ' N \X±, A 2 ) converges weakly to M(Ai,A 2 ). 
Consequently, for any continuous bounded function 0(A) in II 2 , 

lim / 0(A)rfM^(Ai,A 2 ) = / 0(A)dM(Ai,A 2 ). 

N->ocJ U 2 J U 2 

M(Xi, A 2 ) is called the regression spectral measure of X t in accordance with 
the time series literature (Taniguchi et al. 2008). 

(e) is somewhat restrictive, but (e) holds for a random process with the 
separable covariance functions such as the product of 1-dimensional Matern 
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covariance functions (Yajima and Matsuda 2008). (h) is often called axial 
symmetry. The spectral density functions of a separable covariance and 
isotropic one satisfy (h). 

Next we state the asymptotic covariance matrix of three estimators, that 
is Pglse, Ppbe and $lse- 

Define 

D N 2 = diagiWxxW, . . . , \\x p \\) 

where x { = (^(i,i),i, • • • , x (i,N),i, x (2,i),i, ■ ■ ■ , x {2,N),i, ■ ■ ■ , x (N,i),i, ■ ■ ■ , x (N,N),i) 
and Hxjll = (E*T=i E^=i(^i,* 2 ),i) 2 ) 1/2 (i = 1, • • • ,P)- 

Theorem 1 (Yajima and Matsuda (2008)). Under (a)-(g), 



lim D N 2E 



Af->oo 



0> 



GLSE 



GLSE 



D N 2 



(27r) 2 y n2 /(A 1 ,A 2 



■dM(X 1 , A 2 



Theorem 2 (Yajima and Matsuda (2008)). Under (a)-(d) and (f), 



lim D^E 

iV-j-oo 



D 



N 2 



(Plse - P)(Plse - (3) 
= (2tt) 2 R(0, 0)- 1 [ /(Ai, X 2 )dM(X 1 , X 2 )R(0, 0)" 1 . 

Moreover, Yajima and Matsuda (2008) gives the following necessary and 
sufficient conditions where the asymptotic covariance matrix of the LSE is 
identical with that of the GLSE. 



Theorem 3 (Yajima and Matsuda (2008)). Under (a)-(g), the LSE is 
asymptotically efficient relative to the GLSE (BLUE) if and only if M(Ai, A2) 
increases at not more than p values of A (e [0, n] 2 ) and the sum of the ranks 
of the increases in M(X 1 , A 2 ) is p. 

Theorem 1-3 can be regarded as an extension of Grenander and Rosen- 
blatt (1957) for d — 1 to spatial processes. The following theorem is our 
main result. The proof is given in Appendix B. 
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Theorem 4. Under (a)-(d) and (f)-(h), 



lim D^E 



= (2tt) S 



{PPBE - f3)(PPBE - P) 

dM(Xi, A 2 ) 



D 



N 2 



n 2 <7(Ai, A 2 ) 



V/ /(Ai, 



/(Ai,A^ 
A 2 ) 



dM(Ai,A 2 ; 



x 



/ -7TVr^(Ai,A 2 ) 
in 2 0(Ai, A 2 ) 



Note that conditions of Theorems 1-4 are somewhat different. In particu- 
lar, unlike Theorems 2 and 4, Theorems 1 and 3 do not hold due to (e) if the 
random process has the isotropic Matern covariance function. The separable 
covariance functions can satisfy all the conditions of Theorems 1-4. 

Amemiya (1973) and Engle (1974) investigated the asymptotic properties 
of (3pBE f° r d = 1 case and Theorem 4 is an extension of their results to the 
spatial case. 



4 Computational experiments 

We conduct Monte Carlo simulation by using MATLAB. The finite accuracy 
and the asymptotic variance of the PBE are compared with those of the LSE 
through computational simulations especially in the case where the LSE is 
not asymptotically efficient. We consider the linear regression models with 
one regressor 

for 1 < ti,t 2 < N and (3 = 2. For the regressor X( tl ,t 2 ), the polynomial trend 
X(t!,t 2 ) — tit2, the harmonic trend X{t lt t 2 ) = cos((7r/2)ii) cos((7r/2)t 2 ) and the 
polynomial plus harmonic trend X( tl)t2 ) = 1 + cos((7r/2)ti) cos((7r/2)t 2 ) (see 
Toyooka 1985 for d — 1 case) are considered. These regressors satisfy (a)- 
(d) and (g). The jumps of M(Ai,A 2 ) are 1 at (Ai,A 2 ) = (0,0) and 1/4 at 
(Ai,A 2 ) = (7r/2,7r/2),(-7r/2,7r/2),(7r/2,-7r/2),(-7r/2,-7r/2) for the poly- 
nomial trend and the harmonic trend respectively. Therefore, for these two 
regressors the LSE is asymptotically efficient from Theorem 3 if (e) and (f) 
hold. Additionally, it follows from a routine calculation that the asymptotic 
variance of the PBE is identical with that of the LSE in just these cases. As 
a result, the PBE is also asymptotically efficient under the above situation. 
In the third regressor, the jumps of M(Ai, A 2 ) are 4/5 at (Ai, A 2 ) = (0, 0) and 
1/20 at (Ai, A 2 ) = (tt/2, tt/2), (-tt/2, tt/2), (tt/2, -tt/2), (-tt/2, -tt/2). Since 
the condition of Theorem 3 is not satisfied, the LSE is not asymptotically 
efficient. 
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For the true covariance function j € (h), we consider the following six mod- 
els. The first and second models are 



where 



le (h) = C(\\h\\) \\h\\ = yjhl + hl 



and K u (-) is the modified Bessel function of the second kind of order v. This 
is the isotropic Matern covariance function and we put v = 2, p = 3, a 2 = 1 
and z/=l,p = 3,<r 2 = l respectively. The third model is the product of 
1-dimensional Matern covariance functions with v = 2, p = 3, a 2 = 1 and 
v = l,p = 3, o" 2 = 1 respectively. The fourth one is the product of 1- 
dimensional Matern covariance function with u = l,p = 3,a 2 = l and the 
autocovariance function of AR(2) 



a 2 f 2 f 2 

(x) 



gl-\x\ gl-\x\ 



where ^ = (2/3)(l + V3i),^,2 = (2/3)(l — v^^) and a 2 chosen such as 
c*(0) = 1. It is inappropriate to approximate this autocovariance function 
c*(x) by that of AR(1) because the successive negative correlation exists. 
The fifth one is the product of the autocovariance function of AR(2) in the 
fourth one and that of AR(1) 

2 

c**(x) = j^<P lxl , 

where = 0.5 and a 2 ^ chosen such as c**(0) = 1. The sixth one is the 
product of two autocovariance functions of AR(1) with = 0.9 and the scale 
parameter such as 7 e (0) = 1. The covariance functions except the first and 
second ones satisfy (e) (Yajima and Matsuda 2008) and all the models satisfy 
(f) and (h). 

Next we explain the separable approximation like the Yule- Walker esti- 
mator in this simulation. Three types, that is AR(l)xAR(l), AR(l)xAR(2) 
and AR(2)xAR(2), are adopted. The first approximation is expressed by 

70i,/i 2 ) = 71(^1)72(^2) 

2 "> 
Ci ,1ft, 1 cr. 



fell °2 Jh 2 \ 
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Each parameter is estimated by 

7(1,0) 2 7(0,1 



where 



and 



<t>i 



7(0,0)' 



and af 9 

7(0,0) 



7(0,0) 



7' (O,O;0!,0 2 ) 



m 

S(h) 



1 



JV(h) 



(ij)eS(h) 

{(i,j)|ft = ti - t,, ti, t,- G P N }, N(h) = #{S(h)}, 

j AT AT 

%i,t 2 ) = J/(ti,t a ) - $LSEX( tl ,t2h i = Yl Yl 



,*2) 



t 1 = l t 2 =l 



7 0,0; 



(l-«) (l- 



(see Cressie 1993). <5" 2 2 is the estimator of of of and is not necessary for 
the calculation of $pbe- It is used in the estimation of the spectral density 
function g(Ai,A 2 ) as we discuss later. For the second one, the covariance 
function used in the separable approximation and the estimator are 



l(hi,h 2 ) = 71(^1)72(^2 



1 - 



C2^1 



t 1— l^al 
S2 



- a+ \/a 2 + 46 , - 

6 = -j and £ 2 

—lb 



£i 2 -l £1-1 
a/o 2 + 4b 



-lb 



where 



p 2 (0) = 1, p 2 (l) 



P 2 (0) p 2 (-l) 

p 2 (i) p 2 (0) 

7(0,1 



7(0,0) 



p 2 (l) 

p 2 (2) 



p 2 (-l) and p 2 (2) 



7(0,2) 
7(0,0)' 



0! and ofof are estimated by the similar way in the case of the first approx- 
imation model. As before, the parameters are estimated in AR(2)xAR(2) 
case. 



10 



Table 1: Simulation results in the polynomial trend case. 





case 


Matern [y 


= 2.0) 


Matern 


(y = 1.0) 


(is = 2.0) 


x(u = 1.0) 


approximation 


N 


LSE 


PBE 


LSE 


PBE 


LSE 


PBE 


AR(1) x AR(1) 


20 


1.328 


1.130 


1.214 


1.027 


1.187 


1.154 




60 


1.108 


1.123 


1.137 


1.022 


1.066 


1.109 


AR(1) x AR(2) 


20 


_ 


1.073 


_ 


1.032 


_ 


1.134 




60 


_ 


1.057 


_ 


1.021 


_ 


1.082 


AR(2) x AR(2) 


20 




1.038 




1.031 




1.026 




60 




1.018 




1.025 




0.999 




case 


{y = 1.0)x 


AR(2) 


AR(1) 


xAR(2) 


AR(1) 


xAR(l) 


approximation 


N 


LSE 


PBE 


LSE 


PBE 


LSE 


PBE 


AR(1) x AR(1) 


20 


1.145 


1.199 


1.066 


1.179 


1.929 


1.085 




60 


1.073 


1.102 


1.053 


1.091 


1.379 


0.999 


AR(1) x AR(2) 


20 




1.024 




1.006 




1.156 




60 




1.009 




1.0 




1.013 


AR(2) x AR(2) 


20 




1.007 




1.016 




1.182 




60 




1.003 




1.001 




1.019 



We generate the error terms {^(t 1 ,t 2 )}i<t 1 ,t 2 <N by the multivariate normal 
distribution with the mean and the covariance matrix S. (3pbe is calculated 
by using the parameters estimated from €( tl j 2 ) — U(ti,t 2 ) ~ Plsex^m)- The 
sample size is A^ 2 = 20 2 or 60 2 with 1000 replications. From them, we 
calculate the empirical ratios of the variances of the LSE and the PBE to 
that of the GLSE. 

Tables 1 and 2 summarize the results in the case of the first and second 
regressors. As A" increases, the finite efficiency seems to go to 1 in all the 
cases. The PBE shows good performance even in the case of the isotropic 
Matern class. If AR(l)xAR(2) and AR(l)xAR(l) are the true covariance 
function, the PBE by the corresponding approximation is superior to those 
in other cases. Moreover, although the asymptotic variances of the LSE and 
the PBE are equal to each other, the PBE outperforms the LSE in many 
cases of A^ = 60. That is because the PBE includes the information of the 
approximated spatial correlation structure. Throughout these simulations, 
it seems that as the order of the autoregressive process in the separable 
approximation is larger, the efficiency tends to be often better. 

Table 3 shows the result in the polynomial plus harmonic trend case. (•) 
means the theoretical ratios of the asymptotic variances of the LSE and the 
PBE in Theorems 2 and 4 to that of the GLSE in Theorem 1. g(\\, A 2 ) with 
the average of values of parameters estimated in each iteration in A" = 60 
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Table 2: 


Simu 


ation results 


in the harmonic trend case. 








case 


Matern [y 


= 2.0) 


Matern 


{y = 1.0) 


{y = 2.0) 


x(z/ = 1.0) 


approximation 


N 


LSE 


PBE 


LSE 


PBE 


LSE 


PBE 


AR(1) x AR(1) 


20 


1.178 


1.033 


1.040 


1.025 


1.483 


1.049 




60 


1.064 


1.008 


1.032 


1.003 


1.151 


1.015 


AR(1) x AR(2) 


20 


- 


1.035 


- 


1.038 


- 


1.044 




60 


- 


1.012 




1.001 


- 


1.017 


AR(2) x AR(2) 


20 




1.045 




1.048 




1.024 




60 




1.012 




1.003 




1.003 




case 


{y = 1.0)x 


AR(2) 


AR(1) 


xAR(2) 


AR(1) 


xAR(l) 


approximation 


N 


LSE 


PBE 


LSE 


PBE 


LSE 


PBE 


AR(1) x AR(1) 


20 


1.343 


1.106 


1.203 


1.078 


1.770 


1.010 




60 


1.119 


1.033 


1.075 


1.032 


1.271 


1.0 


AR(1) x AR(2) 


20 




1.017 




1.006 




1.016 




60 




1.003 




1.0 




1.0 


AR(2) x AR(2) 


20 




1.010 




1.011 




1.012 




60 




1.0 




1.0 




1.003 



Table 3: Simulation results in the polynomial plus harmonic trend case. 





case 


Matern (y 


= 2.0) 


Matern 


(y = 1.0) 


(y = 2.0)) 


<{v = 1.0) 


approximation 


N 


LSE 


PBE 


LSE 


PBE 


LSE 


PBE 


AR(1) x AR(1) 


20 


37.347 


1.024 


18.060 


1.069 


54.357 


1.051 




60 


42.525 


1.036 


20.083 


1.022 


65.797 


1.028 






(44.431)* 


(1.012)* 


(21.069)* 


(1.019)* 


(70.077) 


(1.004) 


AR(1) x AR(2) 


20 




1.066 




1.070 




1.076 




60 




1.036 




1.029 




1.016 








(1.01)* 




(1.017)* 




(1.003) 


AR(2) x AR(2) 


20 




1.085 




1.099 




1.105 




60 




1.011 




1.002 




1.006 








(1.006)* 




(1.01)* 




(1.002) 




case 


{y = 1.0) > 


AR(2) 


AR(1) 


xAR(2) 


AR(1)> 


<AR(1) 


approximation 


N 


LSE 


PBE 


LSE 


PBE 


LSE 


PBE 


AR(1) x AR(1) 


20 


3.263 


1.260 


1.661 


1.364 


1691 


16.18 




60 


3.652 


1.137 


1.660 


1.327 


3732 


1.02 






(3.575) 


(1.147) 


(1.622) 


(1.283) 


(5242) 


(1.0) 


AR(1) x AR(2) 


20 




1.023 




1.049 




32.515 




60 




0.996 




1.005 




1.022 








(1.002) 




(1.0) 




(1.0) 


AR(2) x AR(2) 


20 




1.068 




1.069 




68.844 




60 




1.002 




1.016 




1.03 








(1.001) 




(1.002) 




(1.001) 
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Table 4: The time required for the calculation of each estimator for three 
approximations (seconds). 



LSE 


GLSE PBE(AR(l)xAR(l)) 


PBE(AR(l)xAR(2)) 


PBE(AR(2)xAR(2)) 


estimation procedure (sec.) 


0.0096 


0.0086 


0.013 


calculation time of estimators (sec.) 4.43 * 10~ 4 


384.793 2.368 


2.404 


2.399 



case is used when we calculate the asymptotic variance of Ppbe in Theorem 
4. Since the isotropic Matern case does not satisfy (e) and the asymptotic 
variance of the GLSE in this case cannot be calculated by Theorem 1, (■)* 
means the ratios of the asymptotic variances of the LSE and the PBE in 
Theorems 2 and 4 to that of the GLSE calculated formally in Theorem 1. 
Note that the spectral density function of the Matern covariance function 
must be calculated in consideration of the aliasing phenomenon (Fuentes 
2005) because the observations are in the integer lattice. Since the LSE is 
not asymptotically efficient in the third regressor case, the performance of 
the LSE is too bad. However, the efficiency of the PBE is near 1 in both the 
empirical and theoretical ratios. As A" increases, the finite efficiency seems 
to go to the value of (•) and (•)*. This supports the theoretical results and 
suggests that (e) may be weakened for Theorem 1. Other properties are 
similar to the preceding simulations. 

Next we see one realization of the spectral density function g(Ai, A 2 ) by 
using the average of values of parameters estimated in each iteration for 
the polynomial plus harmonic trend case. Figure 1 shows the true spectral 
density function / and the approximated spectral density functions g on 
[0,7r] 2 because / and g are axial symmetry. In the case of the separable 
approximation AR(2) x AR(2), except the true model of AR(1) x AR(1), 
the fitting of the spectral density function g for / is very good at the points 
(Ai, A 2 ) where M(Ai, A 2 ) has jumps. It is consistent with the result of Table 
3 which shows the good performance of the separable approximation AR(2) 
x AR(2). 

Finally we examine the computational time of (3glse, Plse and (3pbe- 
We put A" = 100 and adopt the Matern covariance function with v — 2, p — 
3, ex 2 = 1 for the polynomial plus harmonic trend. All computations are 
carried out on Linux powered 3.33GHz Xeon processor with 8 Gbytes RAM. 
From Table 4, the LSE is computationally efficient, but Table 3 shows too 
bad performance in this case. The computation time of the PBE including 
the estimation procedure is faster than that of the GLSE and we observe the 
good performance of the PBE in Table 3. 



13 



5 Conclusion and future studies 



Our study proposes the PBE by using the separable approximation of the 
covariance function on a lattice data as an alternative estimator of the GLSE 
which is infeasible due to the computational burden and the unknown covari- 
ance structure in practice. We derive the asymptotic covariance matrix of 
the PBE and see the effect of the misspecification of the covariance function 
for the the GLSE. The PBE by the separable approximation works well in 
many simulations even if the true covariance function is the isotropic one. 
In particular, when the LSE is not asymptotically efficient, the PBE shows 
a better performance than that of the LSE. Moreover, the PBE reduces the 
computational time substantially because of the separable structure. 

The theoretical comparison of the asymptotic covariance matrix of the 
PBE and that of the LSE is a future work. In addition, the extension of the 
true process to the strongly dependent random field should be considered. 

Appendix A : Technical Lemmas 



To derive the asymptotic covariance matrix of (3 pee, we shall prove the 
two lemmas. 



Lemma A. 1. Let (7r^' m2 , 7r™2' m2 > • • • , ^uv™ 2 > 
be the solution of 



mi,m 2 \' 

)■■■■> n N,l J • • • J "iVJV I 



mi,rri2 mi,rri2 7711,7712 _ 

■ • • 1 7T \r 1 , . . . , 7T 



/ 7711, 7712 \ 

/ ^1,1 \ 



\ 7771,7772 I 

\"JV,JV / 



/ j(mi,m 2 ) 
r y(m 1 ,m 2 + 1) 

7(7711, m 2 + N — 1) 

7(77li + l,77l 2 ) 



\ 



\7(mi + iV - 1, m 2 + N - I) J 
For fixed mi > —N + 1 and m 2 > — iV + 1, 

N N 



\ 11 \ 11 I 7711,7712 I , 

suPZ^ZJ 7 ^' l <0 °- 



7 = 1 j=l 



Proof. flA.ll) is expressed by 



Si 



/ 7771 ,7772 \ 



7771,7772 . 



I 71 (mi) ^ 
^yi(m x + N-l)J 



mi, m 2 G Z. (A.l) 



/ 7 2 (m 2 ) ^ 
\7 2 (m2 + N-1)J 
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Then, from the property of the Kronecker product (see Horn and Johnson 
1991), 







/ 


/ 7i(™i) ^ 







( 

^2 


( l2(m 2 ) 


\ 


\ 


mi ,m 2 1 
\ n N,N / 




V 


\7i( m i + N — 1) J 


J 




V 


\y2(m 2 + N - 




J 


Now let (vr^,vr^ fc ,.. 


■ > ^tv fc)' be the solution of 











\ ir mk I 
for k = 1,2. We will prove 

TV 

sup 22 K\fc I < 00 ( m k > ~N + 1) 



\lk{m k + N 



N 



(A.2) 



i=i 



is 1 



by the mathematical induction. In —N + 1 < m k < (k = 1,2), it" 

^*,(— m fc +i) from the uniqueness of (7r^, 7r^, . . . , vr^ fc fc )' where 5, 
if % = —rrik + 1, otherwise 0. Therefore, (jA.20 is clear in this case. When 
m k = 1, 

/7 fe (l)\ 



If is AR(Pfc) with the autocovariance function 7^, that is = <pi,kZt-i + 
h <pp k:k zt-p k + r) t where % is a white noise, 7rj fe = <p j)k (j = 1, . . . , P fc ) and 

71 j k = (j > -Pfc)- Therefore, supjy Yl%=i l^fcl < 00 • Next assume that (IA.2[) 
holds for mk = m. Consider m k = m + 1. Then, from a routine calculation, 



sp( Z i,..., ZJ v)^/v+m+i = ^=1^™ where Pgp( zl ,..., ZN ) is the projection to 

the closed subspace sp(^i, . . • , Zn) spanned by (zi, . . . , z N ). Moreover, from 
the similar calculation, 

Psp(zi,...,ZN)ZN+m+l ^sp(2i,...,zjv) (P > sp(2i,...,2jv+i)^A r +m+l) 



sp(z 1 ,...,z N ) 



E 

1=1 



■n lk z N+2 -i 



N+l 



Psp(z 1 ,...,z N )^l t k z N+l + 2j ^l,k Z N+2-l 

1=2 

I N \ AT+1 



.1=1 



1=2 
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Therefore, vr™ +1 = vr™ 7rJ fc + 7rft lfJfc (/ = 1, . . . , iV). In this case, 

TV JV TV 

su p E l^ff" 1 ! - sup E KfelKfcl + su p E l^fcl < °°- 

dS2D holds for m fc = m + 1. Note that k™ 1 ' 1 "* = O™ 2 for 1 < i,j < N. 
Then 

N N N N 

EEw ma i = Ei<i 1 iEwi 

i=l j=l i=l j=l 

is bounded for N. □ 

Lemma A. 2. Suppose that h(Xi, X2) is a continuous function on [—it, it] 2 
and h(\i, X2) = h(—Xi,X 2 ) = h(Xi,—X 2 ) = h(— Ai, — A 2 ). T/tere, for any 
sufficiently small e > 0, there exist faz(Ai, A2) and hu(Xi, A2) suc/i £/ia£ 

^l(Ai,A 2 ) = a,fc 1 ,fc 2 e * ( ' fclAl+A:2A2 \ a fci,fc 2 = a -fci,-fe25 

fcl=-Xl k 2 =-K 2 
K 1 K 2 

hu(X h X 2 )= E ^i^^ 1 " 1 ^ 2 " 2 ^ 6fci,* 2 = b-Jfcx.-Jfca, 

k 1 =-K 1 k 2 =-K 2 

h L (X u X 2 ) < h(X u X 2 ) < hu(X u X 2 ) 

and 

hu(Xi, A 2 ) - fc L (Ai, A 2 ) < e, (Ai, A 2 ) e [-tt, tt] 2 . 

Proo/. We set A^fcfr, jfc) = J2\ M < ni £|; 2 |<n 2 e«*»+*»> = A^) Aj^) 
where = E|j,|<„, e * ii2/i (* = 1> 2 )- A^G/i) is the Dirichlet kernel, 

2m + 1 if y< = 0. 



DnXVi) 



Additionally, consider 

SnijuJ 1 ^!-, A 2 ) = (^j e ji,j 2 ) e ji,j 2 

|j'l|<m 1.72 1 <"2 



^2 r r h{X 1 ,X 2 ) E E e iil( ^- Al) e^(^- A2 )dAxdA 2 



(27T) 
1 




|jl|<m |j2|<n2 

/i(Ai, A 2 )D ni (A' 1 - Ai)-D n2 (A 2 - A 2 )dAi(iA 2 , 
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where e jld2 = e ^+^ and (h,e juj2 ) = j\ fl n h(\ u X^e-^^dX.dX^ {2uf . 

By defining /i( A i, A2) = h(X\ +2im X) A2 + 27m 2 ) (n l5 n 2 e Z), it can be rewrit- 
ten as 

S nun2 h(\[, A 2 ) = -—x / / /i(Ai - Ai, A' 2 - A 2 )-Dni(Ai)-Dn 2 (A 2 )rfAidA2. 



Moreover, 

«i — 1 712—1 

nin 2 



^ Til 1 77-2 1 /»7f /»7T 

— Si 1>i2 /i(Ai, A 2 ) = / / /i(Ai-Ai, A 2 -A 2 )i^ni(Ai)i ; r„ 2 (A2)dAidA2, 



i 1= i 2 =0 

where K n .{y^) is defined by 



27m 



— E^(lfc) (< = 1,2) 



j*=o 



and is called the Fejer kernel. Finally, it follows from the argument extended 
from Theorem 2.11.1 of Brockwell and Davis (1991; page 69) and the axial 
symmetry of /i(Ai, A 2 ) that for any e > 0, 



m— 1 U2—1 



h(\ 1 ,\ 2 ) Y] y~] S ilti2 h(X u A 2 ) 



nin 2 f- , 

11=0 12=0 



< e, 



uniformly on [— 7r,7r] 2 for sufficiently large ni and ra 2 . In that case, 

ni— lU2 — 1 711—1 712—1 

nin 2 



^ 711 — 1 712 — 1 711—1 712—1 

~ y Sj^kjx^ x 2 ) = y ^^h^ 



i 1= 12=0 



Ji=-(ni-l) J2=-(ri2-l) 



where Cj lj2 = C-j u j 2 = Cj u -j 2 = C-j u -j 2 from the axial symmetry of fo(Ai, A2). 
Therefore, c^^'s are real and 

ni — l 712 — 1 rii — 1 n 2 — 1 

c jiJ 2 ejij2 = Yl Yl c iij2( co s(ji^i+J2A2)+«sin(j 1 Ai+j 2 A2)). 

ji = -(7ii-l) j2=-(7i2-l) Ji=-(ni-l) i2=-("2-l) 

Now it follows that 

ni—l ri 2 — 1 

Y sinO'iAi + j 2 A 2 ) = 0. 

ji=-(m-i) j2=-("2-i) 
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Therefore, E^m-i £ ba 

|<7i2— i c ii>i2 e jij2 i s rea l- For any e > 0, 

ni — 1 712 — 1 m— 1 "2 — i 

^ ^ c jlJ2 e ilJ2 -- < /i(Ai, A 2 ) < ^ ^ c ji,j2 e h,h+2 

Ji=-(ni-l) j2=-(n2-l) j'i=— (ni— 1) J2=— (ri2 — 1) 

for sufficiently large ni and n 2 . Therefore, setting A x = ri\ — 1,K 2 = n 2 — 
1, a felifc2 = 6 fcl ,fc 2 = c fclifc2 except a , = c , - e/2, 6 ,o = c , + e/2, the proof is 
completed. 

□ 

Appendix B : Proof of Theorem 4 

We prove Theorem 4 by the arguments extended from Anderson (1971) 
and Yajima (1994). 

Proof of Theorem 4 



D N 2E 



0PBE-P)0PBE-P)' 



D N 2 



= [D-lX't^XD-^j 1 [D-lx't'^t^XD-^j [D-lx't^XD'l) 
From Theorem 1 the first and third terms converge to 

[D-lX't^XD-S)- 1 -* (JL f^-J-^iM^M^ 
as N — > oo. Here Put 

/(A 1 ,A 2 ) = /( ^ 1 '^ ) x g (A 1 ,A 2 ) 
= h(X 1: X 2 ) x #(Ai, A 2 ). 

Firstly consider the case of 

Li L 2 

MAi,A 2 )= E b h Me t{hXM \ b luh = b^ l2 . 

h=—L\ l 2 =—L 2 

In this case, 

e i((n-ki)\iHh-k 2 )x 2 ) f ^ x ^ dXidX2 



W J — 7T 



= E E 6 'i.fa T re i ^ 1 -* 1+Il ) Al+ ^-*» +fa ) A3 ^(A 1 ,A 2 )dA 1 dA 2 . 

, _ r 7 _ r J— 7T J — 7T 



= — Li l'2 = — L2 
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Therefore, 

where 
S(/i,/ 2 ) = 

7 (Z 1 + jV-l,Z 2 ) 



S= E E ^(Zi.fe), 

h=—L\ h=—L2 



i{hM) 
7 (Z 1 + JV-l,Z 2 -l) 



7 (Z 1 + 1-7V,Z 2 ) 
7 (Z 1 + 1-7V,Z 2 + 1) 

7(Zi, fe) 



Vt(Zi + iV - 1, Z 2 + JV - 1) 7 (Zi + TV - 1, Z 2 + TV - 2) • • • 7 (/i, Z 2 + TV - 1) 

• • • 7 (Zi + 1 — iV, / 2 + 2 — iV) 7 (Zi + 1 - N, Z 2 + 1 - N) \ 
■■■ 1 (l 1 + l-N,l 2 + 3-N) 1 (l 1 + l-N,l 2 + 2-N) 



7 (l u l 2 + 2-N) 



7 (Z!,Z 2 + 1 -N) 



i(k,k + i) i{hM) J 

Consider the case of l\ > and l 2 > for sufficiently large N. In that case, 
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the (to, n)th element (D^X't^tih, ^t^XD]^) is 



1 



1/2 ) ^(l.O.W ' ' ' ' ^(JV,-) 



(affi^O, 0)^(0, 0)^ 

x ( Z(h+l,l 2 +l),n-> • • • J s (Ji+l,JV),n3 0, . . . , 0, . . . , X(7v,« 2 +l),n ; • • • , •E(N,N),n > 0, . . . , 0, 

0,...,0) + (a;(i 5 .) im , • • • , ac^.^jjjS 1 



x 



AT AT AT AT 

n n V^V^ h+l-JV,l h+l-N,l 2 

'---' '2^Z^ 7r M ^(A r +l-i,A r +l-i),r l ,---,2^Z^ 7r iJ 
i=l j=l i=l j=l 

N N N N 

X X(AT + i_ iiA r + i_j) in , . . . , 0, . . . , 0, 2J ^2 n i'j X (N+l-i,N+l-j),n, • • • , /J ^i'j 2 

i=l j=l i=l j=l 

N N 

x a;(jv+i-i,jv+i-i),nj ^(Ar+i-i.Ar+i-i),™, • • • , 

i=i j=i 

7r 'j ; ' 2a; (JV+l-i,iV+l-j),n I > (B.l) 
i=l i=l / J 

where ac(y ), m = 0(;,i), m , • • • , a;(»,jv),m)' (« = 1, • • • > N ) an d {ttIj 1 '" 12 } is defined 
in Lemma [A. II It follows from the property of Kronecker product (see Horn 
and Johnson 1991), the exact form of the inverse matrix of AR(Pj) (i = 1, 2), 
(a)-(d) and a routine calculation that the first term of flB.lj) converges to 



1 f Q-'i-ih^i+hM) 

dM mn (X 1 ,X 2 ) 



(2tt)2 J n2 g l (X 1 )g 2 (X 2 ) 
as A" — > 00. Next 

1 



I the second term of (lB.ip | < 



a^ N \0,0)a^ N \W 1 ' 2 



X {(^(1, 0.mJ • • • ' X (N,-),m)^ 1 ( £C (l,-),m3 • • • J X (N,-),m) } 

V V \ / N N 

°> ■ ■ ■ ' / J / ! 7r'j /2 a;(v+i-i,v+i-i),n 1 S" 1 1 0, . . ■ ,7^7^ ^/ 2 ^(Ar+i-i,Ar+i-j),n 

t=i j=i / V t=i j=i 
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N N 



1/2 



< 



i,j),m 



aS V) (0,0)air ) (0,0)) 1/2 V U U 



{/ h / N N 

~ Xn2 [ Yl E EEC'^^+Hw-An 
\mi=li+l-)Vm 2 =l \i=l j=l 

h k / N N \ 2 \ ~\ 1 / 2 

+ EE^'^^+w-j)," 

mi=l 7Tl2=Z2+l— iV \ z=l j = l 

1 /2 

,r / , JVmax i<'..fa<Af4.,w.- \ Irlfi-i 



mi=l \i=l 

2 



+ E Era + E E EEKr 2 i) f. (b.2) 

m 2 =l \j=l / mi=lm 2 =l \i=l j=l / J 

where is the greatest eigenvalue of E _1 and (vr^, . . . , 7r™ fe fe ) (/c = 1, 2) is 

defined in the proof of Lemma lA. 1[ Since it is proved that Ajy2 is bounded 
from the property of Kronecker product (see Horn and Johnson 1991) and 
Proposition 4.5.3 of Brockwell and Davis (1991), (IB.2I) converges to as 
N — > oo from (g), Lemma [A. II and its proof. Therefore, 

V J m ,n (2ir) 2 J U 2 g 1 {\i)g 2 {\2) 

as A" — > oo. In a similar way, we can show other cases of (l\ > 0,l 2 < 0), 
(h < 0, h > 0) and (h < 0, l 2 < 0). Then 

V / m,n V / m,n 

l\=—L\ l 2 =—L2 

1 r /(A 1; A 2 ) 

^(27py n 2?(A^) dMm " (Ai ' A2) 

as iV -)■ 00. Finally consider the case of general /i(Ai,A2). From (f), (h) 
and Lemma [A. 21 for any sufficiently small e > 0, there exist h^Xi, X 2 ) and 
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^[/(Ai,A 2 ) such that 

K x K 2 

^l(Ai,A 2 )= ^ ^ afei,A:2 el( ' fclAl+fc2A2 ' > ; a fci,fc 2 — a -fci,-fc 2 , 
fci=-ifi k 2 =-K 2 
K! K 2 

^£/(Ai,A 2 ) = ^ ^ frfci,fc 2 e l ^' lAl+fe2A2 \ b klt k 2 = b-k lt -k 2 , 

fei=-Ki k 2 =~K 2 

/il(Ai,A 2 ) < MAi,A 2 ) < /ic/(Ai,A 2 ) 

and 

^c/(Ai, A 2 ) - /i L (Ai, A 2 ) < e, (Ai, A 2 ) e [-tt, tt] 2 . 
Then for any 7 G K p , 

7 ' j d^x'e- 1 e l s- 1 a;£)^7 < V^-U'E-^E-^D-h < yD-U'E-^E-^D-h, 

where and E^ are the covariance matrices with the spectral density func- 
tions 

^l(Ai, A 2 )y(Ai, A 2 ) and hu(\i, A 2 )(?(Ai, A 2 ) respectively. From the above dis- 
cussion, 

^'7^ / ^MTdMiXuXfrK luni'D-lx't-^XD~h 
(2tt) 2 J n2 c/(Ai, A 2 ) tv^oo 

< ImT V^X'S^EE^X^h < l / M^rfM(A 1) A 2 ) 7 . 

iv->oo (2tt) 2 J n2 p(Ai, A 2 ) 

Since e is arbitrary, as TV — > 00 

7^-U'E-EE-^-h -+ 7'^ £ ^^(A, A 2)7 

(Svr) 2 V r(Ai, A 2 ) 
Since this holds for every vector 7, as iV — >• 00 

The proof is completed. 

□ 
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Figure 1: Approximated spectral density functions on [0,7r] 2 . The true 
covariance function corresponds to the isotropic Matern [y = 2.0), the 
isotropic Matern {y = 1.0) and [y = 2.0) x(u = 1.0) for the left, middle 
and right column in the top row respectively. Similarly, it corresponds to {y 
= 1.0)xAR(2), AR(l)xAR(2) and AR(l)xAR(l) for the left, middle and 
right column in the bottom row respectively. 
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